\section*{The seed length dilemma: long seeds or short seeds} \label{dilemma}

%Although using different data structures and mapping algorithms, the high level
%mapping strategies of both hash-based mappers and suffix-array based mappers are
%very similar. 
%Mappers from both categories index the reference genome into a
%key-value lookup table, as shown in Figure~\ref{fig:Lookup}, where the keys are
%sequences of DNA and the values are lists of appearance locations of each key in
%the reference genome.  mapping, mappers from both categories extract one to
%multiple sub-sequences from the read, which are {\bf seeds}, and use the seeds
%to query the lookup table for locations. These locations are considered as
%potential mapping candidates of the read and are verified through edit-distance
%calculation.

%STOP TALKING ABOUT BWT-FM FROM THIS POINT ON!!!!!!!!!!!!!

An important parameter for hash-based mappers is the seed length. Most
hash-based mappers use seeds at a length between 10 to 14 bases, with some
exceptions which use longer seeds up to 20 bases causing higher memory usage. 
In this section, we propose the metric
``mapping cost'' that depends on the cheap vs. expensive status of seeds (Figure~\ref{fig:Lookup})
to estimate the amount of computation different seed lengths
lead to and analyze the effects of seed lengths on the speed, memory-efficiency
and error tolerance of the mapper.

\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{figures/Lookup_B.pdf}
\caption{A typical lookup table used in hash-based
mappers, which is a permutation array. Logically, the keys of the table are
seeds and the values of the table
are locations. In this table, most seeds have fewer than 100 locations, which
are defined as ``cheap seeds'' while a few seeds have more than 100 locations
which are defined as ``expensive seeds''. In reality, this is implemented as a
permutation array, stores a
pointer to the location list of all permutations. Permutations that have no location
are defined as ``empty seeds'' which store a NULL pointer.}
\label{fig:Lookup}
\end{figure}

\subsection*{Mapping cost and the effects on the mapping speed of the mapper}
\label{performance}

Seed extension, or verification, through alignment  is
the major calculation in hash-based mappers which occupies more than 90\% of
the execution time through our profiling experiments. 
 The shorter the seed is, the more frequently
will it appear in the reference genome thus the more candidate locations it will return. 
As previous works point out\cite{Xin2013}, most of the verification
computation are unnecessary since most locations of the short
seeds are false locations which do not provide correct mappings of the entire
long read. 

In addition, the number of locations of the seeds are extremely unbalanced:
 while most of the seeds have fewer than a hundred locations,
 some seeds may be seen in thousands of locations. We call such seeds as
 {\bf expensive seeds} (solid shaded entries in Figure~\ref{fig:Lookup}). 

 In order to illustrate how expensive seeds hurt the performance of the mapper
 and to evaluate how longer seeds may alleviate this problem, we propose a new
 metric called the {\bf mapping cost}. The mapping cost is defined as the
 average number of locations of a seed selected randomly from the genome. In
 other words, the mapping cost is the estimation of the number of locations
 $E[Loc(S)]$ of a random seed $S$, which is calculated as the sum of products of
 the number of locations of a seed $S_{i}$ times the probability of selecting
 this seed $S_{i}$ through a random draw in the reference genome,  as shown in
 Equation~\ref{eq:map_cost_1}:

\begin{equation} \label{eq:map_cost_1} E[Loc(S)] = \sum_{i}^{N} P(S_{i}) \times
Loc(S_{i}) \end{equation}

The probability of picking the seed $S_{i}$ from a random locus in the
genome is calculated as the number of locations of the seed $S_{i}$ over the total
number of locations of the entire genome, which is the length of the genome defined as
 a constant $C$: $\sum_{i}^{N} Loc(S_{i}) = C$. As a result,
we can rewrite Equation~\ref{eq:map_cost_1} as Equation~\ref{eq:map_cost_2}:

\begin{equation} \label{eq:map_cost_2} E[Loc(S)] = \sum_{i}^{N} P(S_{i}) \times
Loc(S_{i}) = \sum_{i}^{N} \frac{Loc(S_{i})}{C} \times Loc(S_{i}) = \frac{1}{C}
\times \sum_{i}^{N} Loc(S_{i})^{2} \end{equation}

Our goal is to minimize  $E[Loc(S)]$ for a genome of length
$C$. We can formulate our objective function as a classical mathematical problem  known as the Inequality of
Arithmetic Mean and Root-Mean Square (a direct consequence of the Cauchy-Schwarz
Inequality). The inequality says $\sqrt{\frac{x_{1}^{2} + x_{2}^{2} +
\cdot\cdot\cdot + x_{N}^{2}}{N}} \geq \frac{x_{1} + x_{2} + \cdot\cdot\cdot +
x_{N}}{N}$ with equality if and only if $x_{1} = x_{2} = \cdot\cdot\cdot =
x_{N}$. In this case, if all seeds have the same number of locations, in other
words if the locations are distributed evenly into all seeds, we can reach the
minimum of the estimate $E[Loc(S)] = \frac{1}{C} \times N \times
(\frac{C}{N})^{2} =\frac{C}{N}$, which also decreases with a larger $N$. This
implies that we prefer a more balanced location distribution among all seeds and
a larger number of seeds.

\begin{figure}[h]
	\center
	\begin{tabular}{c}
		\subfloat[]{
		\begin{minipage}[c][]{0.8\textwidth}
			\centering
			\includegraphics[width=\textwidth]{figures/Distri_locationsVSseedLen_full_B.pdf}
		\end{minipage}
		} \\
		\subfloat[]{
		\begin{minipage}[c][]{0.8\textwidth}
			\centering
			\includegraphics[width=\textwidth]{figures/Distri_locationsVSseedLen_notFull_B.pdf}
		\end{minipage}
		}
	\end{tabular}
	
	\caption{Figure (a) shows the distribution of seeds in different number of
	locations with different lengths of the human chromosome 1. Figure (b)
	shows the same distribution after taking out the seeds having 0-50
	locations to show the reduction of expensive seeds.}

    \label{fig:seedcompo}
\end{figure}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.7\textwidth]{figures/MappingCostVsSeedLen_B.pdf}
	\caption{Mapping cost with different seed lengths of the human chromosome 1.}
	\label{fig:mapcost}
\end{figure}

%\begin{figure}[h]
%	\subfigure{\includegraphics[width=0.5\textwidth]{figures/Distri_locationsVSseedLen_full_B.pdf}}
%	\subfigure{\includegraphics[width=0.5\textwidth]{figures/Distri_locationsVSseedLen_notFull_B.pdf}}
%	
%	\caption{Distribution of seeds in different number of locations with
%	different lengths of the human chromosome 1. The second plot took out the
%	seeds having 0-50 locations to show the reduction in expensive seeds.}
%    \label{fig:seedcompo}
%\end{figure}
%
%\begin{figure}[h]
%	\centering
%	\includegraphics[width=0.7\textwidth]{figures/MappingCostVsSeedLen_B.pdf}
%	\caption{Mapping cost with different seed lengths of the human chromosome 1.}
%	\label{fig:mapcost}
%\end{figure}

A more balanced location distribution can be accomplished through increasing
the seed length. When the seeds are longer, not only are there more seeds in
total, which is a larger $N$, but also 
the locations are more evenly
distributed among the seeds---all seeds will all together have fewer locations,
generating fewer verification calculations during mapping, as illustrated in
Figure~\ref{fig:seedcompo}. As a result, as the seed length increases, the
mapping cost decreases as shown in Figure~\ref{fig:mapcost}, which improves the mapping speed.

\subsection*{The effects on memory-efficiency of the mapper}
\label{data structure}

As Figure~\ref{fig:seedcompo} also shows, as the seed length increases, the
{\it potential} number of seeds also increases. In fact, the potential number
of seeds at 20 bases is 13.7x of the number of seeds at 12 bases. As a result,
longer seeds require more memory space to store more data in the lookup table,

Theoretically, the size of the lookup table of length 20 should be at most
13.7x of the size of length 12 since the number of seeds of length 20 is 13.7x
of the number of seeds of length 12. However, to perform hash table lookups in
O(1) time, most mappers tend to index all possible permutations of seeds
(Figure~\ref{fig:Lookup}), which exponentially increases the memory
requirements. Among all permutations, many of them never appeared in the
reference genome hence have no locations. In the remainder of the
paper, we call the seeds with zero locations as \textbf{``empty seeds''}.

Empty seeds reduce the memory-efficiency of the mapper since they occupy empty
space in the array. As the seed length increases, a larger portion of the
permutations becomes empty seeds, indicating reduced
memory-efficiency of the mapper, as shown inFigure~\ref{fig:emptyseed}.

% in reality, it is substantially higher than 13.7x, which can
%be nearly 100x or even more. This is because hash based mappers stores not only
%the seeds actually appeared in the reference (i.e., seeds with at least one location) but
%the entire permutations of the bases, starting from AAAAAAAAAAAA to
%TTTTTTTTTTTT, as shown in 
%to enable fast accessing. 
%Note that 
%the upper bound for the {\it real} number of seeds is equal to the genome length. 


%If the lookup table stores only the seeds appeared in the reference, then each
%lookup would would need to search the lookup table for the matching key. Given
%the large number of seeds, this search involves non-trivial computation with
%multiple memory accesses. On the contrary, with permutation arrays, a seed
%lookup becomes a simple index hashing calculation which substitutes each base
%in the seed with a number between 0 to 3. After the index is calculated, mapper
%access the array with the index, and the array returns the pointer leading the
%locations of the querying seed.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.7\textwidth]{figures/EmptySeedPercentage_B.pdf}
	\caption{Percentage of empty seeds in the permutation arrays with different
	seed lengths of the human chromosome 1.}
	\label{fig:emptyseed}
\end{figure}

%The downside of permutation arrays is low memory-efficiency. The permutations
%that never appear in the reference genome, which we call them \textbf{empty
%seeds} (line shaded entries in Figure~\ref{fig:Lookup}), still occupy empty
%slots in the array with a NULL pointer. At 12 bases, most permutations are not
%empty seeds. However as the seed length increases, the percentage of empty
%seeds in the permutation array also increases rapidly. At 20 bases, there are
%in total $4^{20}=1099511627776$ permutations, which is even larger than the
%total number of locations of the entire human genome (3.2 billion). As a
%result, at 20 bases, most permutations in the lookup table are empty seeds, as
%shown in Figure~\ref{fig:emptyseed}, drastically reducing the memory-efficiency
%of the mapper.

\subsection*{The effects on error-tolerance of the mapper}
\label{tolerance}


%During lookup, most mappers assume no error in the seeds (except, the {\it spaced seeds}). 
%However, if a read does
%contain minor errors then the mapper can only map this read correctly when the
%mapper can find at least one intact seed which contains no error.

According to the Pigeonhole Principle, $e$ errors can destroy at most
$e$ non-overlapping seeds, therefore we can guarantee to find an intact seed if the read is
divided into more than $e+1$ non-overlapping seeds. 
When the seeds are short, a read can be divided into many non-overlapping
seeds, which yields high error-tolerance. 
 As the seed length increases, the number of seeds to which a read can be
partitioned into decreases, reducing the error-tolerance.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.7\textwidth]{figures/MaxErrorToleranceVsSeedLen_B.pdf}
	\caption{The maximal number of tolerable errors for a 80 bases read with different
	seed lengths.}
	\label{fig:errortolerance}
\end{figure}

Figure~\ref{fig:errortolerance} show that for a 80-base-long read, the number of tolerated
errors decreases with longer seeds.

